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SOME EXPERIMENTAL RESULTS CONCERNING THE 
ERROR PROPAGATION IN RUNGE-KUTTA TYPE 
INTEGRATION FORMULAS 

INTRODUCTION 

1. In two earlier reports [ i] , [2] , the author derived RUNGE-KUTTA 
formulas up to the eighth order. Each of these formulas, in fact, 
represented a pair of RUNGE-KUTTA formulas. By adding one or two 
more evaluations of the differential equation and changing the weight 
factors, the nth-order formula (n ^ 8) was extended to a (n+1) st-order 
formula . 

2. In these earlier reports, the (n+1) st-order formula was used as a step- 
size control for the nth -order formula since the (n+ 1) st-order formula 
covers correctly the leading term of the local truncation error of the 
nth -order formula. 

3. In this report we shall deal with the global error propagation of our 
RUNGE-KUTTA formulas.. The problem will be approached in two 
different ways. In Section I, we will present the more conventional 
approach using the integrated differential equation for the error propaga- 
tion. In Section II, two-sided ( or bilateral) RUNGE-KUTTA formulas 
are derived. Such two-sided RUNGE-KUTTA formulas are convenient 
tor the investigation of the error propagation of extensive systems of 
differential equations, since no partial derivatives of the differential 
equations are required for this type of formulas. 

4. For both approaches the knowledge of the leading term of the local 
truncation error is essential. In the first approach, the local truncation 
error enters directly the equation for the error propagation. For the 
second approach, our RUNGE-KUTTA formulas can be easily converted 
into two-sided formulas by making use of the leading term of the local 
truncation error of our formulas. 

5. As described in Section I and Section n, both approaches for the error 
propagation can serve to obtain realistic upper and lower bounds for the 
error, along with the integration of the differential equations of the 
problem. In general, the true error will lie somewhere between these 
bounds. However, certain conditions have to be satisfied to guarantee 



that the true error is always located between these bounds . If such 
conditions can be formulated, they probably would be of little practical 
help because they would be too complicated and would involve unknown 
partial derivatives. This report will not be concerned with the formula- 
tion of such conditions. Section III, however, will present some nontrivial 
examples to show that our procedures are capable of yielding reasonably 
close error bounds for the solutions of such problems . 


SECTION I. THE ERROR EQUATIONS, BASED ON ONE 
INTEGRATION PROCEDURE PER INTEGRATION STEP 


6. For simplicity, let us consider a system of only two first-order 
differential equations: 


— =f(t,x,y) 

-f =g ( t,*,y> 


( 1 ) 


The numerical integration of ( 1) is performed by a RUNGE-KUTTA 
formula. For the first equation ( i) , such a RUNGE-KUTTA formula 
would read 


f 0 = f(to,Xo,y 0 ) \ 

fj = f(to + c^h, Xq + h/3 10 f 0 , y 0 + h/3 i0 g 0 ) / 

f 2 = f[to + a 2 h, Xq + h(/3 20 fo + ^ 21 fi) , y 0 + h(/3 20 g 0 + ftigi)] > ( 2) 


and 


x = Xq + h(c 0 f 0 + Cift + c 2 f 2 + . . . ) (3) 

The coefficients ai,a! 2 , . . . ; • • • ’> aQ d c 0 ,c 1 ,c 2 , . . * 

are known numerical constants of the RUNGE-KUTTA formula, and 
h(= dt) stands for the integration stepsize. 


2 



A corresponding RUNGE-KUTTA formula holds for the second differential 
equation ( 1) . 


7. We now assume that the values Xo and y 0 at the beginning t = t 0 of 
our integration step are affected by errors ( e x ^ 0 and and we 

like to study the propagation of these errors through the current integra- 
tion step. 


Including these errors in (2) and (3) , we obtain 

V f [VV( £ x)o' V( £ y)o] 

h " f to + a i h ' V ( e x)o + W- y 0 + ( £ y)o + ' S ioV > ] 

*2 = '['0 + “ 2 h ’ X 0 + ( £ x)0 + (Vo + Vl) h ' 

y O + ( £ y)o + (Vo + ^21 ? l) h ] 


(4) 


x + e = + (e \ + h /c f + c .f . + c + . . . ^ 

x 0 \ x/0 \ 0 0 11 22 ) 

If we expend ( 4) in TAYLOR series and carry only linear terms in 
and ( e y) q » ^he following expressions result: 


(5) 


1. We disregard errors in t 0 since such errors can be avoided by a proper 
selection of the time step. For a binary electronic computer, time steps 
that are a (positive or negative) power of 2 should be used. 
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The three subscripts vOO(v = 0, 1,2, ... ) in (6) are supposed to 
indicate that the expression in question is to be taken at 

( = 0) . For example stands for 


i oe ca 

(D 


*o + % h ' x o' y o 

taken at t 0 + a^h; XojYo* 


Introducing (6) into (5) yields an expression for the error e at the end 

of the integration step. If, in this expression, we restrict ourselves to 
linear terms in h, we may drop in (6) all terms that are multiplied with 


h and may replace f’g^yQO by — and 



000 


(-f)voo by (-t)ooo since 


(M) „ (X\ + (-iL.) 

\dx/v00 \dx/000 \dxdtj 


/df\ = /_df\ / 9 2 f \ 

\0y/ rOO ~\ay/ 000 \9y9t/ 


a h + 
000 v 


a h + 
000 V 


(7) 
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Because of the equation of condition 


c 0 + Ci + c 2 + . . . = 1 


( 8 ) 


for the RUNGE-KUTTA coefficients, the expression for 
corresponding expression for e then read 



e and the 
x 


(9) 


In (9) , we have replaced the triple subscript 000 of the partial deriva- 
tives by just one subscript 0. 


Equation (9) represents the propagation of the error 



0 


through the current integration step. To obtain the total error, we still 
have to add to the right-hand sides of (9) the local truncation and local 
round-off error committed during the current integration step. 


Denoting these errors by T , T and R , R , respectively, we 

x y x y 

obtain instead of (9) in vector form 



In an obvious way, equation ( 10) can be extended to systems of more 
than two first-order differential equations. 


8. To handle these well-known equations ( 10) for the error propagation, the 


partial derivatives , -g , -g . g 


have to be computed along with 


the integration of the differential equations. Furthermore, certain 
assumptions have to be made for the local truncation errors T , 1 
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and the local round-off errors R^, R^. In the following, such assump- 
tions for these errors will be discussed. 

9. As for the truncation error, we restrict ourselves to the leading term 

of this error. In our earlier reports [ 1] , [2] , we already determined 

the leading truncation error term for RUNGE-KUTTA formulas up to 

the eighth order. This leading term of the local truncation error was 

used in these earlier reports for the stepsize control only. We introduce 

this term now also in (10) for T and T . 

x y 

10. An accurate determination of the round-off errors R and R in ( 10) 

x y 

is almost impossible. We therefore resort to a somewhat crude but 
easily obtainable approximation. In (3), the quantities x 0 and 
h(c 0 f 0 + c-jf^ + c 2 f 2 + . . . ) are in general of different order of magnitude, 
the latter expression being a small increment of Xg. Therefore, when 
performing their addition in (3) , the electronic computer has to shift a 
certain part of the smaller number h(c 0 f 0 + cjfj + c 2 f 2 + . . . ) out of its 
range. The shifted -out part of the smaller number is lost for the 
computation. We now consider this shifted-out part as an approximation 
for the round-off error R in ( 10) . Naturally, the functions 

X 

f 0 , f lt f 2 , ... in (2) are also affected by round-off errors . However, 
when these functions enter equation ( 3) , their round-off errors are 
multiplied by the small quantity h and may therefore be neglected, 
compared with the round-off error resulting from the addition of x 0 
and h(c 0 f 0 + Cjf* + c 2 f 2 + . . . ) . 

In the same way, an approximate value for the round-off error R^ can 
be obtained from the formula for y that corresponds to (3) . 


11. Introducing the approximate values T , T and R , R , the partial 
... . . x y x y 


derivatives (-§)„, (f*) 0 , (f *) 0 and the total errors (e^, 

^ q at the beginning of the current integration step into the equations 
( 10) , the total errors e^, at the end of the current integration step 
are easily obtained by a few multiplications and additions. 
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We evaluated equations (10) three times: 



/ ("A /*. 

/ € x \ / 9x 9y \ 

\ e ^ / \_9g 9g / 

\ ^ / \3x 3y/0 






(ID 


(v = 0, 1,2). 


For v = 1 we took proper care of the round-off errors R , R , thereby 

( 1) ( 1) x y 

obtaining approximations e , e for the true errors. For 

x y 

v = 0 (disregarding the round-off errors) and for v = 2 (doubling the 

round-off errors) we established "bounds" e ^ , e ^ and e \ 

( 2 ) x y x 

e for the true errors . 

y 

(v) (v) 

Compensating for these errors 


for v = 1 approximations x^ ^ , y^ ^ for the true values x^,, y, 


, we obtained from x, y 

For 
( 2 ) 


T* 
and x' 


v = 0 and y = 2, the compensation yields bounds x^ , y^ 

( 2 ) 

y ’ for x^, and y^. The values x,y were the results of the RUNGE- 

KUTTA integration (2) , (3) of our problem without any error compensa- 
tion. 


Figure 1 illustrates the procedure for the computation of the error spread 
in x for the first two integration steps. 

12. As already stated in the Introduction, we cannot always expect strict 
bounds using such an extremely simplified approximation procedure for 
the round-off errors. However, the examples of Section III show that by 
using in (11) the correct values of the partial derivatives and proper 
approximations for the local truncation errors, realistic and reasonably 
close error bounds are obtainable by our procedure for the round-off 
errors . 
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In some problems, it might occur that our "bounds" are too tight. In 

such cases, the values x^ , x^ , and x^ (or the corresponding 
values of y) would become practically equal, but not necessarily equal 
to the true value, x^,. Better bounds might then be obtained by replacing 

in ( 11) v = 0,1,2 by v= -1,1,3 or similar values to effect a wider 
spread of the bounds. 


13. 


If the partial derivatives in ( 10) assume large values during the numerical 

integration of the problem ( 1) , the error can propagate heavily, even if 

the local truncation errors T and T are kept very small by a 

x y 

sufficiently small integration stepsize. To slow down suchlan undesirably 
large error spread, we introduce an additional test for the stepsize by 


requiring that the products 



h 


and 



should not exceed a pre-given value. 


If they exceed that 


value, the stepsize h will be halved sufficiently often until all products 
stay below the pre-given value. 


By a proper choice of this pre-given test value, we can regulate the error 
spread to a certain extent. However, one should keep in mind that a very 
small test value might lead to very small integration stepsizes . Thereby, 
a heavy buildup of the round-off errors might occur, resulting in 
unrealistically large error bounds. 


SFCTION II. THE ERROR EQUATIONS, BASED ON TWO 
INTEGRATION PROCEDURES PER INTEGRATION STEP 


14. Although the method of Section I, if properly applied, will lead to rather 
accurate values for the errors and reasonably close error bounds, the 
method has the disadvantage of requiring the computation of the partial 


derivatives 



These partial derivatives might turn 


out to be lengthy and cumbersome expressions . Furthermore, for 
extensive systems of differential equations there will be a large number 
of such partial derivatives ( n 2 partial derivatives for a system of n 
differential equations) . 
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15. To avoid these inconveniences in the study of the error propagation, one 
might rather resort to methods which do not involve partial derivatives. 
Such methods can be established if for each step the numerical integration 
is performed twice by means of two different RUNGE-KUTTA formulas. 
These two RUNGE-KUTTA formulas are chosen in such a way that their 
leading terms of the local truncation error are equal but of opposite sign. 

Obviously, the arithmetic mean of these two RUNGE-KUTTA formulas 
represents a RUNGE-KUTTA formula of order n+i, if the two original 
RUNGE-KUTTA formulas are of order n. Considering this arithmetic 
mean formula as an approximation of the true values, the difference in 
the values of the original two formulas can be regarded as an indicator 
for the spread of the truncation error. 

In the literature, two such RUNGE-KUTTA formulas are called two-sided 
or bilateral RUNGE-KUTTA formulas. 

Since such two-sided RUNGE-KUTTA formulas do not consider the round- 
off errors, these must be kept negligible. This means, the error propaga- 
tion obtained from two-sided RUNGE-KUTTA formulas will be reliable 
only as long as the truncation errors are dominant compared with the 
round-off errors. 

Because of the dominance of the truncation errors, the error bounds of 
two-sided RUNGE-KUTTA formulas will not be as close as the bounds 
of the procedure described in Section I. This is not surprising, since the 
truncation errors can be practically eliminated in the procedure of 
Section I. 

16. Two-sided RUNGE-KUTTA formulas of low order have been known a 
long time. For example, the modified EULER-CAUCHY formula 

x= Xo + “• h {f (t 0 ,x 0 ) + f [t 0 + h, x 0 + hf (t 0 ,Xo)]} , (12) 

which represents a RUNGE-KUTTA formula of the second order, can be 
considered as the arithmetic mean of the two first-order RUNGE-KUTTA 
formulas: 


9 



= Xq + hf(t 0 ,x 0 ) ) 

x^ 2 ^ = Xq + hf [t 0 + h, Xq + hf (t 0 ,Xo)] | 

TAYLOR-expansion of the right-hand sides of ( 13) leads to: 


( 13 ) 



Since the last terms on the right-hand sides of ( 14) represent the 
leading terms of the local truncation error of formulas ( 13) , these 
formulas are indeed two-sided RUNGE-KUTTA formulas of the first 
order. 


In the case of formulas ( 12) and ( 13) , the conditions that ( 13) always 
yields strict bounds for the solution ( 12) have been established by 
S. GORN and R. MOORE [3] . These conditions are also quoted in 
another paper by S. GORN ( [4] , p. 76 ) . They involve the first- and 
second-order partial derivatives f , f , f and f . 

t tt tx XX 

More recently, two-sided RUNGE-KUTTA formulas of the second- and 
third -order were published by A.D, GORBUNOV and YU. A. SHAKHOV 
[5], [6]. 

17. Two-sided RUNGE-KUTTA formulas, up to the eighth order, can be 
obtained easily from the RUNGE -KUTTA formulas of our reports [ 1] , 

[2] . We have only to recall that the formulas of these reports represent 
pairs of RUNGE-KUTTA formulas of the order n and n+1 ( 1 <n g 8) . 
The formula of the order n+ 1 is obtained from the formula of the order 
n by adding one or two more evaluations of the differential equations and 
changing the weight factors c of the nth-order formula to the weight 

K 

factors c of the (n+i)st-order formula. 

K 
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Knowing two such formulas of the nth- and the (n+i)st -order, one can 
immediately derive from these formulas two two-sided RUNGE-KUTTA 
formulas of the nth -order by putting 


c 


(i) 


K 


2c - c 
K K 


c <2 > = 

K 


C 

K 


(15) 


However, it is not even necessary to compute the coefficients ( 15) , since 


we can operate our RUNGE-KUTTA formulas of [ 1] , [2] as two-sided 
RUNGE-KUTTA formulas in the following way: Starting from the initial 

values ^t Q ,x 0 ^ , . . . ^ , we compute one step by our nth -order RUNGE 


KUTTA formula applying the leading term TE of the local truncation 
error for the stepsize control. From the result, we subtract twice the 

term TE, thereby obtaining the values x^ ^ , ... for one of our two- 
sided RUNGE-KUTTA formulas. Next, we start from the initial values 


^tg,Xg^ , . . . ^ and compute one step with the same stepsize as for 

^ , thereby obtaining the values x^ for our second two-sided 
RUNGE-KUTTA formula. 


In Figure 2 the procedure is illustrated for the first two integration steps 
The final values x are then approximated by 



and the propagated errors by 



(17) 


Corresponding formulas hold for y, e and for any further dependent 
variables . ^ 

18. For the convenience of the reader, we have listed as Appendix to this 
report our RUNGE-KUTTA coefficients for RK1(2), RK2(3), . . . , 
RK7(8) from [2] , [ 1] , including the leading term TE of the local 
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truncation error. However, formulas RK1(2) and RK2(3) in [2] 
have been changed to formulas more suitable for the study of error 
propagation. For the somewhat lengthy table of coefficients for RK8(9) 
we refer to our paper [ 1 ] . 


19 . 


While the local truncation error can be checked by the term TE of the 


Appendix, we have control of 
integration step by checking 


the error propagation during the current 


e - 


( £ y) 


etc . , since 


'x ( £ x)ol ’ 

these expressions represent the error growth during the current integra- 
tion step (the suffix 0 indicating the values at the beginning of the 
By requiring that the expressions j ^ J , 

etc . remain smaller than a pre-given small quantity and 


current step) 

|£ y-( £ y) 


by reducing the stepsize h until these condition are met, we have, 
similar as in the procedure of Section I, a certain control of the error 
propagation. 


20. To get the error propagation of two-sided RUNGE-KUTTA formulas 

properly started, the local truncation error has to be dominant from the 
beginning of the integration. Therefore, in the beginning, we relax the 
tolerance for the leading term of the local truncation error. After the 
propagated error has grown to a certain magnitude, we continue the com- 
putation by using a sharper (smaller) tolerance. In this way, for the 
examples of Section III, we could keep the error bounds of our two-sided 
RUNGE-KUTTA formulas relatively close (about three to four times as 
large as the bounds for the procedure of Section I) . 


SECTION III. EXAMPLES FOR THE ERROR PROPAGATION 


2i. We apply the methods of Sections I and II to the following examples: 


Problem I: 


dx _ _i 
dt _ 2 


x 

t + i 


2t * y 


dy = J_ . _Z 

dt 2 t + 1 


( 18 ) 


+ 2t • x 


For the initial values 


t 0 = 0, x 0 = i, y 0 = 0 , ( 19) 

the differential equations ( 18) have the following solution in closed 
form: 

X,j, = N/'t + 1 * 

y T = *ft+ 1 • 

Figures 3 and 4 show the shape of the (true) solution (20) . Because of 
the argument t 2 in (20) , the frequency of the oscillations increases rapidly 
with t. The problem clearly requires a continuous stepsize control as 
provided by our RUNGE-KUTTA formulas. We integrated problem ( 18) , 

( 19) from t = 0 to t = 5 . 


cos(t 2 ) 

sin(t 2 ) 


( 20 ) 


Problem n — Restricted Problem of Three Bodies; (see V. SZEBEHELY 
[7] or other textbooks of celestial mechanics): 

In the rotating coordinate system, the problem leads to the well-known 
differential equations: 


dx dy x+jx 

w at [<«+,.)»♦ ✓f* 


- M • 


x - m 1 


[(x-m’) 2 + y 2 ] 3/z 


d 2 y dx , 

d?- = y - 2 •“ 


dt 


[(x+ju) 2 + y 2 f / ' z 


- M 


[(x-/Lt’) 2 + y 2 ] 3//z 


( 21 ) 
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Here, x and y represent the coordinates of a moving particle of 
negligible weight (spaceship) with respect to the barycentric system of 
the two fixed bodies (earth, moon). In this system, the earth has the 
constant coordinates x = - n, y_ = 0, and the moon has the constant 

k -Cj 

coordinates x M = n' = 1 - y M = 0. 


We solved problem (21) for two different sets of initial values: 


a. t 0 - 0, Xo=l. 15, y 0 =O, Xo = 0, y 0 = - 0. 87 I = 


80 


( 22 ) 


and 


b. tb = 0, x 0 = 1. 2, y 0 = 0, xo = 0, y 0 = 




(23) 


Figures 10 and 15 show the shape of the orbit for the initial conditions 
(22) or (23) , respectively. The marks along the orbit represent the 
time scale. The first orbit was computed from t = 0 to t = 6, the 
second orbit from t = 0 to t = 25 . 


Since no solution in closed form for problem (21) , (22) or problem 
(21) , (23) is known, we substituted for the true solution an orbit 
obtained by carrying more than 20 decimal digits on the computer. 2 The 
numerical solution of our problems was obtained by carrying 16 digits 
(IBM-7094) . 

Since, for some parts of its orbit, the spaceship approaches earth or 
moon relatively closely but for other parts is carried far away from these 
bodies, we again need a continuous stepsize control to take proper care 
of all parts of the orbit. 

22. Figures 5 through 9, 11 through 14, and 16 and 17 show our results for 
the propagated errors Ax of Problems I, Ha, and Kb. The errors Ay 
look similar, especially in Problem I, and were therefore omitted. 


2. This part of the computations was carried out by Mr. F. R. Calhoun 
from the Computer Sciences Corporation (CSC) who also developed the 
necessary 30- and 40-digit packages for the IBM-7094 computer. 
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On these graphs we used a few abbreviated expressions that should be 
explained now: 

a. TOL stands for tolerance. The products TOL • | Xq | , 

TOL • |y 0 1 , etc. were used as tolerances for the local 
truncation errors (x 0 ,y 0 , . . . being the initial values in 
x, y, ... for the current step) . The integration stepsize 

h(= dt) was selected in such a manner that the local truncation 
errors | TE | , | TE | , etc . would not exceed the above products. 

x y 

b . TMPDDT stands for test value for the ( absolute) maximum of 
the partial derivatives times dt. This test value serves as a 
control for the error propagation; (see No. 13) . 

c. TEXY stands for test value for the (accumulated) error in x 
or y. As soon as both of these errors in absolute value 
exceed TEXY, we switched over from a relaxed tolerance 
(TOL = 0.1- 10 i2 ) to a sharper tolerance (TOL = 0. 1 • 10 16 ) . 

d. TOLEP stands for tolerance in error propagation. As soon as 

the maximum of le - (e \ I, I e - se \ I , etc. exceeds 
I x V x/0 1 I y y;0| 

the test value TOLEP, the stepsize h(= dt) is halved (if 
necessary repeatedly halved) to prevent a too fast increase in 
error propagation; (see No. 19) . 

In all three problems, we printed our results in equidistant time intervals. 

Li Problems I and Ha we printed for At = , in Problem lib for 

1 8 
At = — . To obtain results for such equidistant time intervals, we had 

Ci 

to adjust our stepsize h(= dt) when overshooting the time interval. In 
the graphs , the results for the printed time intervals were connected by 
straight lines . 

23. For Problem I, the numerical integration was performed by our RUNGE- 
KUTTA formulas RK5(6) and RK7(8). 

Figures 5, 6, and 7 show the results based on the error propagation for- 
mulas of Section I (one integration procedure per integration step) . 
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The curves v = 0, 2 or v = -1, 3 represent the deviations x^ - x^ , 

where x^ ^ stands for the x-value corrected for the local truncation 
and round-off error; (see No. 11 arid No. 12) . The third curve on these 
graphs, which is mostly close to the t-axis represents the deviation 

x^ - x~ ^ between the true solution and x^ ^ 


Naturally, this deviation 


would be zero if the equations for the error propagation would give 
results which are completely correct. 


We see from Figure 5 that for RK5(6) the error bounds increase from 
0 to 4000 • 10 1 = 0.4 • 10 12 . This means that the 13th place behind 
the decimal point could be wrong by up to four units. The real error, 
represented by the third curve in Figure 5, is somewhat smaller than 
our bounds, namely of the order 100 • 10 16 = 0. 01 • 10 2 . We obtain 
safe error bounds in this case. Figure 6 shows the error behavior for 
RK7(8) . Here the error bounds are smaller, since the integration is 
performed in considerably fewer steps. The error bounds do not exceed 


600 


10^ 16 = 0 . 


10 


-13 


However, because of the larger permissible 
^ ^ deviates more from zero than in* Figure 5 . 


stepsize, the curve x^, - x 
It almost reaches the peaks of the error bound curves. 


We recomputed the error behavior for RK7(8) using for the error 
bounds v = -1 and v = 3. The results are plotted in Figure 7. Now the 


curve x T - x^ ^ fits better between the error bounds which now increase 
— 16 — 12 

up to 1200 -10 = 0. 12 • 10 . Again, we obtain safe error bounds in 

this case. 


Figures 8 and 9 show our results based on the error propagation formulas 

of Section II (two integration procedures per integration step) . The 

t (1) , (2) 1 [ (i) (2)1 

curves represent x - x and x - x with x = — [x + x J and 


as third curve, close to the t-axis, x T - x. The error bounds grow up 

to 20 000 • 10 16 = 0.2 • 10 11 for our formula RK5(6) and up to 
6000 • 10 16 =0.6 • 10 for our formula RK7(8). The actual errors 
x T - x are smaller than these error bounds, so that our bounds, again, 


can be considered as safe. 
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Comparing Figures 5, 6, and 7 with Figures 8 and 9, it is very evident 
that the method of Section I leads, as expected, to closer error bounds 
than the method of Section II. All Figures show that our error bounds 
oscillate with approximately the same frequency as the solution itself 
(Fig. 3). 

24. The numerical integration for the orbit of Problem Ila, again, was 
performed by our RUN GE-KUTTA formulas RK5(6) and RK7(8). 
Figures ii to 14 show our results for the error propagation. After the 
explanations of No. 23, our results need hardly any further interpretation. 
It is remarkable how the error bounds grow when the spaceship comes 
close to the earth (t » 1.5 and t ~ 4.875) . Again, our formulas lead 

to safe error bounds, with larger bounds for RK5(6) than for RK7(8) , 
and the error bounds for the method of Section I are closer than those 
for the method of Section II. 

25. Finally, we integrated Problem lib by our RUNGE-KUTTA formula 
RK7( 8) . Because of the longer time range to be covered in this problem 
and because of the more complicated shape of the orbit — the spaceship 
approaches the earth five times — our lower-order formula RK5(6) 
proved not to be very suitable for this problem. Figures 16 and 17 
show the results for our formula RK7(8) . Especially in Figure 17, we 
can easily recognize the increase in the error bounds when the spaceship 
approaches the earth (t « 2, t « 6.5, t ~ 12.5, t » 17.5, and 

t ~ 22.5) . 

26 . Although there is no guarantee that the methods of Sections I and II yield 
strict bounds that bracket the correct solution, we obtained reasonable 
bounds for the problems of this Section. The bounds bracket the correct 
solution almost everywhere, and they are also reasonably close to the 
correct solution. It might require some experience to find the proper 
values for the test functions involved (TMPDDT or TEXY, TOLEP), 
especially in problems ( such as our Problems Ha and lib) in which the 
differential equations contain singularities. In general, it will be 
advisable to use high-order RUNGE-KUTTA formulas, since their 
application will keep the propagated errors small. Because such high- 
order formulas require a relatively small number of integration steps, 
the errors have no real chance to propagate heavily. 


George C. Marshall Space Flight Center 

National Aeronautics and Space Administration 

Marshall Space Flight Center, Alabama 35812, May 26, 1970 
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Figure 1 
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. Error spread for methods using one integration procedure per 
integration step. 



Figure 2. Error spread for methods using two integration procedures per 

integration step. 




Figure 3. Solution x T = n/T+T • cos (t 2 ) for Problem I. 

CO 


I 





Figure 4. Solution y T = 
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Figure 7. Propagated 
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~1 


error Ax for Problem I, using RK7(8) and one integration procedure 
per step (v = -i, 1, 3) . 
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RK5(6) and one integration procedure 
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Figure 16. Propagated error Ax for Problem lib, using RK7(8) and one integration procedure 

per step (v= -1,1,3). 
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